/* Copyright 2001,2002,2003 NAH6
 * All Rights Reserved
 *
 * Parts Copyright DoD, Parts Copyright Starium
 *
 */
#include "main.h"
#include "pctorc.h"

/**************************************************************************
*
* NAME
*       pctorc
*
* FUNCTION
*
*       Convert from lp-polynomial to reflection coefficients.
*
*       BEWARE: This code does not use memory efficiently.
*
* SYNOPSIS
*
*       subroutine PCtoRC(lpc, rc)
*
*   formal
*                       data    I/O
*       name            type    type    function
*       -------------------------------------------------------------------
*       lpc(n+1)        fxpt_16  i      Array of n+1 coefficients
*                                       a(0)+a(1)z**(-1) + a(2)Z**(-2) +
*                                       .... + a(n)z**(-n)
*       rc(n)           fxpt_16 i/o     reflection coefficients (voiced-> +rc1)
*
***************************************************************************
*
* DESCRIPTION
*
*       This routine uses the Levinson recursion to compute reflection
*       coefficients from the LPC coefficients.  The first LPC
*       coefficient is assumed to be 1, and although it is passed
*       to the routine, it is not used in the calculations.
*       Note:  the dimension of the internal array t limits the value
*       of the maximum order.
*
*       CELP's LPC predictor coefficient convention is:
*              p+1         -(i-1)
*       A(z) = SUM   a   z          where a  = +1.0
*              i=1    i                    1
*
*       The sign convention used defines the first reflection coefficient
*       as the normalized first autocorrelation coefficient, which results
*       in positive values of rc(1) for voiced speech.
*
***************************************************************************
*
* CALLED BY
*
*
* CALLS
*
*
*
**************************************************************************/
void PCtoRC(
fxpt_16 lpc[ORDER],             /* 2.13 format */
fxpt_16 rc[ORDER])              /* 2.13 format */
{
        int i, j;
        fxpt_32 num;
        fxpt_16 den;
        fxpt_16 t[ORDER+1], a[ORDER+1];

        for (i = 0; i <= ORDER; i++)
                a[i] = lpc[i];

        for (i = ORDER; i > 1; i--) {
                rc[i-1] = fxpt_negate16(a[i]);
                for (j = 1; j < i; j++) {
                        /*t[i-j] = (a[i-j] + rc[i-1] * a[j]) / (1.0 - rc[i-1] * rc[i-1]);*/
                        num = fxpt_mac32(fxpt_shr32_fast(
                            fxpt_deposit_h(a[i-j]), 2),
                            rc[i-1], a[j]);
                        den = fxpt_sub16(8192,
                            fxpt_mult16_fix(rc[i-1], rc[i-1], 13));
                        t[i-j] = fxpt_extract_l(fxpt_shr32_fast(num/den, 1));
                }
                for (j = 1; j < i; j++)
                        a[j] = t[j];
        }
        rc[0] = fxpt_negate16(a[1]);
}
